Keşifsel Analiz

Amundsen, 1911. DALL-E 3

1. Kategorik Değişkenler

Soru

Chicago Üniversitesi’nde bulunan Ulusal Araştırma Görüş Merkezi (NORC), 1972’den bu yana Genel Sosyal Anketi (GSS) yürütmektedir. Anketin amacı, toplumsal değişimi izlemek ve Amerikan toplumunun artan karmaşıklığını incelemektir. Kişisel tanımlayıcılar arasında yaş, ırk, cinsiyet, medeni durum, gelir, eğitim, siyasi parti ve dini bağlantı gibi özellikler bulunmaktadır. 2000 ile 2014 yılları arasında GSS’ye yanıt verenler arasında medeni durum ve gelir dağılımı nasıldı?

#load tidyverse library
library(tidyverse)

Veri

#import GSS data
data(gss_cat)

#review data structure
str(gss_cat)

## tibble [21,483 × 9] (S3: tbl_df/tbl/data.frame)
##  $ year   : int [1:21483] 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 ...
##  $ marital: Factor w/ 6 levels "No answer","Never married",..: 2 4 5 2 4 6 2 4 6 6 ...
##  $ age    : int [1:21483] 26 48 67 39 25 25 36 44 44 47 ...
##  $ race   : Factor w/ 4 levels "Other","Black",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ rincome: Factor w/ 16 levels "No answer","Don't know",..: 8 8 16 16 16 5 4 9 4 4 ...
##  $ partyid: Factor w/ 10 levels "No answer","Don't know",..: 6 5 7 6 9 10 5 8 9 4 ...
##  $ relig  : Factor w/ 16 levels "No answer","Don't know",..: 15 15 15 6 12 15 5 15 15 15 ...
##  $ denom  : Factor w/ 30 levels "No answer","Don't know",..: 25 23 3 30 30 25 30 15 4 25 ...
##  $ tvhours: int [1:21483] 12 NA 2 4 1 NA 3 NA 0 3 ...

#select two columns and omit rows with missing data
gss_cat2 <- gss_cat %>%
  select(marital,rincome) %>%
  na.omit()

#show levels of marital status
levels(gss_cat2$marital)

## [1] "No answer"     "Never married" "Separated"     "Divorced"     
## [5] "Widowed"       "Married"

#show levels of income
levels(gss_cat2$rincome)

##  [1] "No answer"      "Don't know"     "Refused"        "$25000 or more"
##  [5] "$20000 - 24999" "$15000 - 19999" "$10000 - 14999" "$8000 to 9999" 
##  [9] "$7000 to 7999"  "$6000 to 6999"  "$5000 to 5999"  "$4000 to 4999" 
## [13] "$3000 to 3999"  "$1000 to 2999"  "Lt $1000"       "Not applicable"

#combine and rename income levels
gss_cat3 <- gss_cat2 %>%
  mutate(rincome2=fct_collapse(rincome,
                   missing=c("No answer","Don't know","Refused","Not applicable"),
                   under_5=c("Lt $1000","$1000 to 2999","$3000 to 3999","$4000 to 4999"),
                   from_5to10=c("$5000 to 5999","$6000 to 6999","$7000 to 7999","$8000 to 9999"),
                   from_10to15="$10000 - 14999",
                   from_15to20="$15000 - 19999",
                   from_20to25="$20000 - 24999",
                   over_25="$25000 or more")) %>%
  select(-rincome)

#show new levels of income
levels(gss_cat3$rincome2)

## [1] "missing"     "over_25"     "from_20to25" "from_15to20" "from_10to15"
## [6] "from_5to10"  "under_5"

#reorder income levels
gss_cat4 <- gss_cat3 %>%
  mutate(rincome2=fct_relevel(rincome2,
                   c("missing","under_5","from_5to10","from_10to15","from_15to20",
                     "from_20to25","over_25")))

#show new categories of income
levels(gss_cat4$rincome2)

## [1] "missing"     "under_5"     "from_5to10"  "from_10to15" "from_15to20"
## [6] "from_20to25" "over_25"

Analiz

#summarize variables
summary(gss_cat4)

##           marital             rincome2   
##  No answer    :   17   missing    :8468  
##  Never married: 5416   under_5    :1183  
##  Separated    :  743   from_5to10 : 970  
##  Divorced     : 3383   from_10to15:1168  
##  Widowed      : 1807   from_15to20:1048  
##  Married      :10117   from_20to25:1283  
##                        over_25    :7363

#compute income level percentages
gss_cat4 %>%
  group_by(rincome2) %>%
  summarize(count=n()) %>%
  mutate(prop=count/sum(count)) %>%
  ungroup()

## # A tibble: 7 × 3
##   rincome2    count   prop
##   <fct>       <int>  <dbl>
## 1 missing      8468 0.394 
## 2 under_5      1183 0.0551
## 3 from_5to10    970 0.0452
## 4 from_10to15  1168 0.0544
## 5 from_15to20  1048 0.0488
## 6 from_20to25  1283 0.0597
## 7 over_25      7363 0.343

#plot basic bar chart of marital status
plot(gss_cat4$marital)

#plot ggplot2 bar chart of marital status
ggplot(data=gss_cat4,aes(x=marital)) +
  geom_bar(color="black",fill="sky blue") +
  labs(title="General Social Survey (2000-2014)",
       x="Marital Status",
       y="Count of Respondents",
       caption="Source: University of Chicago (NORC)") +
  theme_bw()

#plot ggplot2 bar chart of income
ggplot(data=gss_cat4,aes(x=rincome2)) +
  geom_bar(aes(y=after_stat(prop),group=1),
           color="black",fill="sky blue") +
  labs(title="General Social Survey (2000-2014)",
       x="Reported Income ($K)",
       y="Proportion of Respondents",
       caption="Source: University of Chicago (NORC)") +
  theme_bw()

Değerlendirme

#plot descending bar chart of marital status
ggplot(data=gss_cat4,aes(x=fct_infreq(marital))) +
  geom_bar(color="black",fill="sky blue") +
  labs(title="General Social Survey (2000-2014)",
       x="Marital Status",
       y="Count of Respondents",
       caption="Source: University of Chicago (NORC)") +
  theme_bw()

#plot misleading bar chart of income
ggplot(data=gss_cat4,aes(x=rincome2)) +
  geom_bar(aes(y=after_stat(prop),group=1),
           color="black",fill="sky blue") +
  labs(title="General Social Survey (2000-2014)",
       x="Reported Income ($K)",
       y="Proportion of Respondents",
       caption="Source: University of Chicago (NORC)") +
  theme_bw() +
  coord_cartesian(ylim = c(0.01,0.99))

Yanıt

Anket yanıtlarında medeni durum ve gelir seviyesini nasıl tanımlamalıyız? 3.4 Şekline bakıldığında, ankete yanıt veren insanların büyük çoğunluğunun evli olduğunu görüyoruz. Evli olmayanların çoğu hiç evlenmemiş ya da boşanmış durumda. Ayrı yaşayan veya dul olan yanıtlayanların sayısı nispeten azdır. 3.3 Şekline göre, yaklaşık %40’lık bir kesim gelir bildirmemişken, yaklaşık %35’lik bir kesim yılda 25 bin dolardan fazla gelir bildirmiştir. Geri kalan %25’lik kesim, 0 ila 25 bin dolar arasındaki beş gelir seviyesine oldukça eşit bir şekilde dağılmıştır.

2. Çaprazlık Tablosu (Contingency Table)

Soru

GSS verilerini kullanarak, daha önce Amerikalılar arasında medeni durumun bireysel dağılımını inceledik. Ancak, medeni durum ve diğer demografik faktörler arasındaki ilişki de ilgimizi çekebilir. 2000 ile 2014 yılları arasında GSS’ye yanıt verenler arasında medeni durum ile siyasi parti arasında bir ilişki var mı?

Veri

#import GSS data
data(gss_cat)

#review marital levels
levels(gss_cat$marital)

## [1] "No answer"     "Never married" "Separated"     "Divorced"     
## [5] "Widowed"       "Married"

#review political levels
levels(gss_cat$partyid)

##  [1] "No answer"          "Don't know"         "Other party"       
##  [4] "Strong republican"  "Not str republican" "Ind,near rep"      
##  [7] "Independent"        "Ind,near dem"       "Not str democrat"  
## [10] "Strong democrat"

#wrangle GSS data
gss_cat2 <- gss_cat %>%
  select(marital,partyid) %>%
  na.omit() %>%
  filter(marital != "No answer",
         !partyid %in% c("No answer","Don't know","Other party")) %>%
  transmute(marital=case_when(
      marital=="Married" ~ "Now",
      marital %in% c("Separated","Divorced","Widowed") ~ "Before",
      marital=="Never married" ~ "Never"),
    political=case_when(
      partyid %in% c("Strong republican","Not str republican") ~ "Republican",
      partyid %in% c("Ind,near rep","Independent","Ind,near dem") ~ "Independent",
      partyid %in% c("Strong democrat","Not str democrat") ~ "Democrat")) %>%
  mutate(marital=as.factor(marital),
         political=as.factor(political))

#summarize data frame
summary(gss_cat2)

##    marital           political   
##  Before:5774   Democrat   :7177  
##  Never :5264   Independent:8403  
##  Now   :9885   Republican :5343

Analiz

#create contingency table
table(gss_cat2$marital,
      gss_cat2$political)

##         
##          Democrat Independent Republican
##   Before     2252        2237       1285
##   Never      1975        2406        883
##   Now        2950        3760       3175

Değerlendirme

#create stacked bar chart
ggplot(data=gss_cat2,aes(x=political,fill=marital)) +
  geom_bar() +
  labs(title="General Social Survey (2000-2014)",
       x="Political Party",
       y="Count",
       fill="Marital Status",
       caption="Source: University of Chicago (NORC)") +
  theme_bw()

#create standardized stacked bar chart
ggplot(data=gss_cat2,aes(x=political,fill=marital)) +
  geom_bar(position="fill") +
  labs(title="General Social Survey (2000-2014)",
       x="Political Party",
       y="Proportion",
       fill="Marital Status",
       caption="Source: University of Chicago (NORC)") +
  theme_bw()

#create standardized stacked bar chart
ggplot(data=gss_cat2,aes(x=marital,fill=political)) +
  geom_bar(position="fill") +
  labs(title="General Social Survey (2000-2014)",
       x="Marital Status",
       y="Proportion",
       fill="Political Party",
       caption="Source: University of Chicago (NORC)") +
  theme_bw()

Yanıt

2000 ile 2014 yılları arasında GSS’ye yanıt verenler arasında, medeni durum ile siyasi parti arasında bir ilişki olduğu görülüyor. Cumhuriyetçiler arasında şu anda evli olanların oranı, diğer iki partiden daha yüksek. Öte yandan, medeni durumun dağılımı Bağımsızlar ve Demokratlar için yaklaşık olarak aynı görünüyor. Ters perspektiften bakıldığında, şu anda evli olanlar, evli olmayanlara (veya hiç olmamış olanlara) göre daha sık Cumhuriyetçi görünüyor. Her iki şekilde de bakıldığında, bu iki faktör ilişkilidir.

Cevabımız, kapsam 2000 ile 2014 yılları arasında sınırlı kaldığı sürece, orijinal soruyla iyi bir uyum içindedir. Ancak, siyasi iklimdeki, sosyal normlardaki ve ekonomik faktörlerdeki değişiklikler, analizimizin sonuçlarını değiştirebilir. Sonuç olarak, tavsiyemiz büyük ölçüde değişebilir. Güncel bir politik kampanyanın yöneticisi için, en güncel rehberliğe erişim kritik öneme sahiptir. Aynı anket sorularına modern yanıtlar sağlandığında, analizimizi kolayca tekrarlayabilir ve yeni içgörüler sunabiliriz.

3. Metin Madenciliği

#load tidyverse and tidytext libraries
library(tidyverse)
library(tidytext)

Soru

Bu vaka çalışmasında, Lewis Carroll (Carroll 1865) tarafından yazılan klasik roman Alice Harikalar Diyarında’yı inceliyoruz. Carroll’ın gerçek adı Charles Lutwidge Dodgson’dı ve Oxford Üniversitesi’nde matematikçi olarak görev yapmanın yanı sıra yazar ve şairdi. Çalışmalarını özellikle geometri, doğrusal cebir ve matematiksel mantık alanlarında yürüttü. Ancak, birçok animasyon ve canlı aksiyon filmine ilham veren bu kitabın yazarı olarak en çok tanınır. Alice Harikalar Diyarında adlı eserde kelime kullanımı ve duygu dağılımı nasıldır?

Veri

#import text file
alice <- read_delim("alice.txt",delim=" ",col_names=FALSE)

## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
##   dat <- vroom(...)
##   problems(dat)

## Rows: 2477 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: " "
## chr (2): X1, X2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

#display text structure
glimpse(alice)

## Rows: 2,477
## Columns: 2
## $ X1 <chr> "CHAPTER", "Down", "Alice", "bank,", "the", "conversations", "“with…
## $ X2 <chr> "I.", "the Rabbit-Hole", "was beginning to get very tired of sittin…

#clean text file
alice2 <- alice %>%
  unite("text",X1:X2,sep=" ",na.rm=TRUE) %>%
  mutate(text=str_squish(text)) %>%
  filter(!str_detect(text,"\\*")) %>%
  transmute(line=1:2468,text)

#display structure of clean data
glimpse(alice2)

## Rows: 2,468
## Columns: 2
## $ line <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19…
## $ text <chr> "CHAPTER I.", "Down the Rabbit-Hole", "Alice was beginning to get…

#tokenize text file by word
alice3 <- alice2 %>%
  unnest_tokens(output=word,input=text,token="words")

#display tokenized structure
glimpse(alice3)

## Rows: 26,687
## Columns: 2
## $ line <int> 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4,…
## $ word <chr> "chapter", "i", "down", "the", "rabbit", "hole", "alice", "was", …

Analiz

#count word frequency
alice3 %>%
  count(word,sort=TRUE) %>%
  head(n=10)

## # A tibble: 10 × 2
##    word      n
##    <chr> <int>
##  1 the    1643
##  2 and     871
##  3 to      729
##  4 a       632
##  5 she     538
##  6 it      527
##  7 of      514
##  8 said    460
##  9 i       393
## 10 alice   386

#load stop words
data(stop_words)

#remove stop words
alice4 <- alice3 %>%
  anti_join(stop_words,by="word")
  
#count frequency without stop words
alice4 %>%
  count(word,sort=TRUE) %>%
  head(n=10)

## # A tibble: 10 × 2
##    word        n
##    <chr>   <int>
##  1 alice     386
##  2 time       71
##  3 queen      68
##  4 king       61
##  5 don’t      60
##  6 it’s       57
##  7 i’m        56
##  8 mock       56
##  9 turtle     56
## 10 gryphon    55

#plot bar chart of word frequency
alice4 %>%
  filter(word!="alice") %>%
  count(word,sort=TRUE) %>%
  head(n=15) %>%
ggplot(aes(x=reorder(word,-n),y=n)) +
  geom_bar(color="black",fill="sky blue",
           stat="identity") +
  labs(title="Alice in Wonderland",
       subtitle="by Lewis Carroll",
       x="Word",
       y="Count",
       caption="Source: Project Gutenberg") +
  theme_bw()

#import sentiment data
sentiments <- get_sentiments("bing")

#join sentiments to words
alice5 <- alice4 %>%
  inner_join(sentiments,by="word")

#summarize sentiments
alice5 %>%
  group_by(sentiment) %>%
  summarize(count=n()) %>%
  ungroup()

## # A tibble: 2 × 2
##   sentiment count
##   <chr>     <int>
## 1 negative    729
## 2 positive    335

#compute cumulative sentiment
alice6 <- alice5 %>%
  mutate(tot_pos=cumsum(sentiment=="positive"),
         tot_neg=cumsum(sentiment=="negative")) %>%
  group_by(line) %>%
  summarize(max_pos=max(tot_pos),
            max_neg=max(tot_neg)) %>%
  ungroup()

#plot cumulative sentiment
ggplot(data=alice6,aes(x=line)) +
  geom_line(aes(y=max_pos),color="blue") +
  geom_line(aes(y=max_neg),color="red") +
  labs(title="Alice in Wonderland",
       subtitle="by Lewis Carroll",
       x="Line Number in Text",
       y="Cumulative Words",
       caption="Source: Project Gutenberg") +
  annotate(geom="text",x=2000,y=675,label="Negative",color="red") +
  annotate(geom="text",x=2000,y=350,label="Positive",color="blue") +
  theme_bw()

Değerlendirme

#load wordcloud library
library(wordcloud)

## Loading required package: RColorBrewer

#create wordcloud of frequencies
alice4 %>%
  filter(word!="alice") %>%
  count(word,sort=TRUE) %>%
  with(wordcloud(words=word,freq=n,
          max.words=100,random.order=FALSE))

#Top 5 negative words
alice5 %>%
  filter(sentiment=="negative") %>%
  group_by(word) %>%
  summarize(count=n()) %>%
  ungroup() %>%
  arrange(desc(count)) %>%
  head(n=5)

## # A tibble: 5 × 2
##   word      count
##   <chr>     <int>
## 1 mock         56
## 2 poor         27
## 3 hastily      16
## 4 mad          15
## 5 anxiously    14

Değerlendirme

Alice hikayesindeki en yaygın olumsuz kelime açık ara ile ‘mock’ (alay etmek) kelimesidir. Bu kelime, yerleşik ‘bing’ sözlüğü tarafından olumsuz duygu olarak atfedilmiştir. Birini alay etmek, onunla kinayeli bir şekilde dalga geçmek anlamına gelir. Sahte veya gerçek dışı olan bir nesne için ‘mock’ (sahte) kelimesi kullanılır. Dolayısıyla, olumsuz duygu atfı mantıklıdır. Ancak, bir alan uzmanı hızla, Alice hikayesindeki bir karakterin Mock Turtle (Sahte Kaplumbağa) olduğunu belirtecektir. ‘mock’ ve ‘turtle’ (kaplumbağa) kelimelerinin 3.9 Şekilde yan yana görünmesi bir tesadüf değildir. ‘Mock’ kelimesi yalnızca Mock Turtle’a atıfta bulunulduğunda kullanılır. Bu bağlamda olumsuz olarak değerlendirilmeli midir? Bu soruya cevap vermek için alan uzmanlarıyla işbirliği gerekebilir.

Mock Turtle örneğinin temel noktası bağlamın önemidir. Analizimiz, bireysel kelimelerin ötesine geçip, duyguyu doğru bir şekilde atamak için diğer kelimeler ve ifadelerle yakınlığını göz önünde bulundurmalıdır. Bununla birlikte, bireysel kelimeler bile bağlama bağlı olarak zıt duygular ifade edebilir. ‘Bad’ (kötü), ‘sick’ (hasta/kötü) ve ‘wicked’ (kötü/harika) gibi kelimeler bazı durumlarda olumsuz, diğerlerinde ise olumlu duygu ifade etmek için kullanılır. Bizim kullanımımıza sunulan çeşitli güçlü metin madenciliği araçları var, ancak dil karmaşık ve inceliklidir. Doğal dil işleme (NLP) alanı geniş olarak bu zorlukları çözmeye çalışır.

Yanıt

Alice Harikalar Diyarında’daki anlamlı kelimelerin sıklığına dayanarak, baş karakter Alice üzerinde anlaşılabilir bir şekilde önemli bir odaklanma bulunmaktadır. Alice’in ötesinde, metin bir dizi ek karaktere atıfta bulunur. Bunlar (azalan sıklıkta) kraliçe, kral, kaplumbağa, grifon, şapka ustası, tavşan ve fareyi içerir. Kitapta görünen önemli temalar zaman, ses ve görme (bakıldı) gibi görünmektedir. Kitap ayrıca sıkça başa atıfta bulunuyor. Hikayeyi tanıyanlar için bu, kraliçenin idam takıntısı nedeniyledir!

Duygu açısından, kitap baştan sona tutarlı bir şekilde daha olumsuzdan pozitife doğru ilerliyor. Bununla birlikte, aşırı olumsuz temanın tutarlı olduğundan emin olmak için duygu analizini başka sözlüklerle (sadece “bing” dışında) tekrar yapmalıyız. Ayrıca, analizi yalnızca pozitif ve negatiften öte duyguları tanımlayan daha ayrıntılı sözlüklerle genişletebiliriz. Diğer genişletmeler, daha güçlü NLP araçları kullanarak çok kelime ve ifade analizini içerebilir. Bu tür takip çalışmaları, “Mock Turtle” (Sahte Kaplumbağa) ve “Off with their heads” (Başlarından edin) gibi bitişik kelimeleri ortaya çıkararak kalan kafa karışıklığını çözebilir.

4. Sayısal Değişkenler

Soru

Major League Baseball (MLB), 1871 yılından bu yana Amerika’da profesyonel beyzbol sporunu düzenlemektedir. Sonuç olarak, beyzbol, mevcut olan en eski spor verilerinden birini sağlar. Gerçekten de, spor analitik endüstrisi, büyük ölçüde Michael Lewis’in Moneyball (Lewis 2004) kitabında anlatılan olaylar nedeniyle popülerlik kazandı. Bu kitap daha sonra büyük bir sinema filmine uyarlandı. Bu tür analizlere kısa bir giriş yapmak gerekirse, şu soruyla karşılaştığımızı varsayalım: 2010 ve 2020 yılları arasında, major lig beyzbolcuları arasında fırsat ve verimlilik dağılımı nasıldı?

Veri

#load Lahman library
library(Lahman)

#review available data sets
data()

#import batting data
data(Batting)

#inspect structure
glimpse(Batting)

## Rows: 112,184
## Columns: 22
## $ playerID <chr> "abercda01", "addybo01", "allisar01", "allisdo01", "ansonca01…
## $ yearID   <int> 1871, 1871, 1871, 1871, 1871, 1871, 1871, 1871, 1871, 1871, 1…
## $ stint    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ teamID   <fct> TRO, RC1, CL1, WS3, RC1, FW1, RC1, BS1, FW1, BS1, CL1, CL1, W…
## $ lgID     <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ G        <int> 1, 25, 29, 27, 25, 12, 1, 31, 1, 18, 22, 1, 10, 3, 20, 29, 1,…
## $ AB       <int> 4, 118, 137, 133, 120, 49, 4, 157, 5, 86, 89, 3, 36, 15, 94, …
## $ R        <int> 0, 30, 28, 28, 29, 9, 0, 66, 1, 13, 18, 0, 6, 7, 24, 26, 0, 0…
## $ H        <int> 0, 32, 40, 44, 39, 11, 1, 63, 1, 13, 27, 0, 7, 6, 33, 32, 0, …
## $ X2B      <int> 0, 6, 4, 10, 11, 2, 0, 10, 1, 2, 1, 0, 0, 0, 9, 3, 0, 0, 1, 0…
## $ X3B      <int> 0, 0, 5, 2, 3, 1, 0, 9, 0, 1, 10, 0, 0, 0, 1, 3, 0, 0, 1, 0, …
## $ HR       <int> 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0…
## $ RBI      <int> 0, 13, 19, 27, 16, 5, 2, 34, 1, 11, 18, 0, 1, 5, 21, 23, 0, 0…
## $ SB       <int> 0, 8, 3, 1, 6, 0, 0, 11, 0, 1, 0, 0, 2, 2, 4, 4, 0, 0, 3, 0, …
## $ CS       <int> 0, 1, 1, 1, 2, 1, 0, 6, 0, 0, 1, 0, 0, 0, 0, 4, 0, 0, 1, 0, 0…
## $ BB       <int> 0, 4, 2, 0, 2, 0, 1, 13, 0, 0, 3, 1, 2, 0, 2, 9, 0, 0, 4, 1, …
## $ SO       <int> 0, 0, 5, 2, 1, 1, 0, 1, 0, 0, 4, 0, 0, 0, 2, 2, 3, 0, 2, 0, 2…
## $ IBB      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ HBP      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ SH       <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ SF       <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ GIDP     <int> 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 1, 2, 0, 0, 0, 0, 3…

#wrangle batting data
Batting2 <- Batting %>%
  filter(yearID>=2010,yearID<=2020) %>%
  group_by(playerID,yearID) %>%
  summarize(tAB=sum(AB),
            tH=sum(H)) %>%
  ungroup() %>%
  filter(tAB>=80) %>%
  mutate(BA=tH/tAB)

## `summarise()` has grouped output by 'playerID'. You can override using the
## `.groups` argument.

#check for impossible batting average
sum(Batting2$BA>1)

## [1] 0

Analiz

#summarize total at-bats
Batting2 %>%
  summarize(count=n(),
            mean=mean(tAB),
            stdev=sd(tAB))

## # A tibble: 1 × 3
##   count  mean stdev
##   <int> <dbl> <dbl>
## 1  4916  326.  168.

#summarize batting average
Batting2 %>%
  summarize(count=n(),
            mean=mean(BA),
            stdev=sd(BA))

## # A tibble: 1 × 3
##   count  mean  stdev
##   <int> <dbl>  <dbl>
## 1  4916 0.250 0.0372

#plot basic histogram of at-bats
hist(Batting2$tAB)

#plot ggplot2 histogram of at-bats
ggplot(data=Batting2,aes(x=tAB)) +
  geom_histogram(color="black",fill="sky blue") +
  labs(title="Major League Baseball (2010-2020)",
       x="Total At-Bats per Season",
       y="Count of Players",
       caption="Source: Sean Lahman (SABR)") +
  theme_bw()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#determine bin number
bins <- sqrt(nrow(Batting2))

#determine range of at-bats
range <- max(Batting2$tAB)-min(Batting2$tAB)

#determine bin widths
range/bins

## [1] 8.614518

#create tailored histogram of at-bats
atbats <- ggplot(data=Batting2,aes(x=tAB)) +
  geom_histogram(binwidth=10,boundary=80,closed="left",
                 color="black",fill="sky blue") +
  labs(title="Major League Baseball (2010-2020)",
       x="Total At-Bats per Season",
       y="Count of Players",
       caption="Source: Sean Lahman (SABR)") +
  theme_bw() +
  scale_x_continuous(limits=c(80,700),breaks=seq(80,700,40)) +
  scale_y_continuous(limits=c(0,150),breaks=seq(0,150,25))

#display histogram
atbats

#create tailored histogram of batting average
batavg <- ggplot(data=Batting2,aes(x=BA)) +
  geom_histogram(aes(y=0.005*..density..),
                 binwidth=0.005,boundary=0,closed="left",
                 color="black",fill="sky blue") +
  labs(title="Major League Baseball (2010-2020)",
       x="Batting Average per Season",
       y="Proportion of Players",
       caption="Source: Sean Lahman (SABR)") +
  theme_bw() +
  scale_x_continuous(limits=c(0.1,0.4),breaks=seq(0.1,0.4,0.05)) +
  scale_y_continuous(limits=c(0,0.07),breaks=seq(0,0.07,0.01))

#display histogram
batavg

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

Değerlendirme

#display at-bats histogram with stats
atbats +
  geom_vline(xintercept=326,color="red",size=1) +
  geom_vline(xintercept=c(158,494),color='red',size=1,linetype="dashed")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#display batting average histogram with stats
batavg +
  geom_vline(xintercept=0.25,color="red",size=1) +
  geom_vline(xintercept=c(0.21,0.29),color="red",size=1,linetype="dashed")

## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

Yanıt

MLB oyuncularının vuruş fırsatını ve etkinliğini nasıl tanımlamalıyız? Fırsat açısından bakıldığında, iki farklı vuruşçu grubu olduğu görünüyor. İlk grup oyuncular, ya sakatlık ya da kötü performans nedeniyle nispeten az sayıda vuruş fırsatına sahip. İkinci bir grup ise sağlıklı kalmayı başarmış ve/veya iyi performans göstererek birçok vuruş fırsatı yakalamış görünüyor. Bu grubun zirve frekansı, 162 maçlık bir sezon üzerinden yayılan yaklaşık 550 vuruş fırsatı etrafında toplanmış. Etkinlik açısından bakıldığında, ortalama vuruş başarısı yüzde 25’tir ve daha düşük veya daha yüksek değerler giderek daha az sıklıkta görülür. Nispeten az sayıda vurucu, sezonu yüzde 20’nin altında veya yüzde 30’un üzerinde bir etkinlikle tamamlar.

5. Ayrık Değerler

Soru

Çocuk Sağlığı ve Gelişimi Araştırmaları (CHDS), 1959’dan beri sağlık ve hastalıkların nesiller arası nasıl geçtiğini araştırmaktadır. Araştırmacıları, yalnızca genetik yatkınlığa değil, aynı zamanda sosyal ve çevresel faktörlere de odaklanmaktadır. Uzun vadeli bir çalışmanın bir parçası olarak, CHDS, bekleyen ebeveynlerin demografik bilgilerini ve sigara içme alışkanlıklarını topladı. Bize şu soru sorulduğunu varsayalım: CHDS projesi kapsamında, hamilelik sırasında sigara içen ve içmeyen annelerin bebeklerinde herhangi aşırı veya beklenmedik doğum ağırlıkları oldu mu?

Veri

#load library
library(mosaicData)

#import gestation data
data(Gestation)

#display structure of data
glimpse(Gestation)

## Rows: 1,236
## Columns: 23
## $ id        <dbl> 15, 20, 58, 61, 72, 100, 102, 129, 142, 148, 164, 171, 175, …
## $ plurality <chr> "single fetus", "single fetus", "single fetus", "single fetu…
## $ outcome   <chr> "live birth", "live birth", "live birth", "live birth", "liv…
## $ date      <date> 1964-11-11, 1965-02-07, 1965-04-25, 1965-02-12, 1964-11-25,…
## $ gestation <dbl> 284, 282, 279, NA, 282, 286, 244, 245, 289, 299, 351, 282, 2…
## $ sex       <chr> "male", "male", "male", "male", "male", "male", "male", "mal…
## $ wt        <dbl> 120, 113, 128, 123, 108, 136, 138, 132, 120, 143, 140, 144, …
## $ parity    <dbl> 1, 2, 1, 2, 1, 4, 4, 2, 3, 3, 2, 4, 3, 5, 3, 4, 3, 3, 2, 3, …
## $ race      <chr> "asian", "white", "white", "white", "white", "white", "black…
## $ age       <dbl> 27, 33, 28, 36, 23, 25, 33, 23, 25, 30, 27, 32, 23, 36, 30, …
## $ ed        <chr> "College graduate", "College graduate", "HS graduate--no oth…
## $ ht        <dbl> 62, 64, 64, 69, 67, 62, 62, 65, 62, 66, 68, 64, 63, 61, 63, …
## $ wt.1      <dbl> 100, 135, 115, 190, 125, 93, 178, 140, 125, 136, 120, 124, 1…
## $ drace     <chr> "asian", "white", "white", "white", "white", "white", "black…
## $ dage      <dbl> 31, 38, 32, 43, 24, 28, 37, 23, 26, 34, 28, 36, 28, 32, 42, …
## $ ded       <chr> "College graduate", "College graduate", "8th -12th grade - d…
## $ dht       <dbl> 65, 70, NA, 68, NA, 64, NA, 71, 70, NA, NA, 74, NA, NA, NA, …
## $ dwt       <dbl> 110, 148, NA, 197, NA, 130, NA, 192, 180, NA, NA, 185, NA, N…
## $ marital   <chr> "married", "married", "married", "married", "married", "marr…
## $ inc       <chr> "2500-5000", "10000-12500", "5000-7500", "20000-22500", "250…
## $ smoke     <chr> "never", "never", "now", "once did, not now", "now", "until …
## $ time      <chr> "never smoked", "never smoked", "still smokes", "2 to 3 year…
## $ number    <chr> "never", "never", "1-4 per day", "20-29 per day", "20-29 per…

#check date max and min
max(Gestation$date)

## [1] "1965-09-10"

min(Gestation$date)

## [1] "1964-09-11"

#wrangle gestation data
Gestation2 <- Gestation %>%
  transmute(weight=wt/16,
            smoke=as.factor(smoke)) %>%
  na.omit()

#display smoke levels
levels(Gestation2$smoke)

## [1] "never"                   "now"                    
## [3] "once did, not now"       "until current pregnancy"

Analiz

#summarize gestation data
summary(Gestation2)

##      weight                           smoke    
##  Min.   : 3.438   never                  :544  
##  1st Qu.: 6.812   now                    :484  
##  Median : 7.500   once did, not now      :103  
##  Mean   : 7.470   until current pregnancy: 95  
##  3rd Qu.: 8.188                                
##  Max.   :11.000

#compute 1st quartile
quantile(Gestation2$weight,probs=0.25)

##    25% 
## 6.8125

#compute 3rd quartile
quantile(Gestation2$weight,probs=0.75)

##    75% 
## 8.1875

#manually compute IQR
8.1875-6.8125

## [1] 1.375

#directly compute IQR
IQR(Gestation2$weight)

## [1] 1.375

#threshold for low outliers
quantile(Gestation2$weight,probs=0.25)-1.5*IQR(Gestation2$weight)

##  25% 
## 4.75

#threshold for high outliers
quantile(Gestation2$weight,probs=0.75)+1.5*IQR(Gestation2$weight)

##   75% 
## 10.25

#filter and count outliers
Gestation2 %>%
  filter(weight<4.75 | weight>10.25) %>%
  nrow()

## [1] 31

#display box plot of birth weight
ggplot(data=Gestation2,aes(x=weight)) +
  geom_boxplot(fill="sky blue") +
  labs(title="Child Health and Development Studies (CHDS)",
       subtitle="Live births between 1964-1965",
       x="Birth Weight (pounds)",
       caption="Source: Public Health Institute") +
  scale_x_continuous(limits=c(3,12),breaks=seq(3,12,0.5)) +
  theme_bw() +
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank())

#display box plot of birth weight by smoking status
ggplot(data=Gestation2,aes(x=weight,y=smoke)) +
  geom_boxplot(fill="sky blue") +
  labs(title="Child Health and Development Studies (CHDS)",
       subtitle="Live births between 1964-1965",
       x="Birth Weight (pounds)",
       y="Mother's Smoking Status",
       caption="Source: Public Health Institute") +
  scale_x_continuous(limits=c(3,12),breaks=seq(3,12,0.5)) +
  theme_bw()

Değerlendirme

#remove all outliers
Gestation2_wo <- Gestation2 %>%
  filter(weight>=4.75, weight<=10.25)

#compute median for both data sets
median(Gestation2$weight)

## [1] 7.5

median(Gestation2_wo$weight)

## [1] 7.5

#compute mean for both data sets
mean(Gestation2$weight)

## [1] 7.469923

mean(Gestation2_wo$weight)

## [1] 7.498431

#compute IQR for both data sets
IQR(Gestation2$weight)

## [1] 1.375

IQR(Gestation2_wo$weight)

## [1] 1.375

#compute stdev for both data sets
sd(Gestation2$weight)

## [1] 1.137723

sd(Gestation2_wo$weight)

## [1] 1.031594

Gestation %>%
  filter(wt/16 <4.5) %>%
  transmute(gestation_months=gestation/30) %>%
  na.omit() %>%
  arrange(gestation_months)

## # A tibble: 12 × 1
##    gestation_months
##               <dbl>
##  1             6.8 
##  2             7.43
##  3             7.6 
##  4             7.73
##  5             7.73
##  6             7.8 
##  7             7.87
##  8             7.9 
##  9             8.17
## 10             8.47
## 11             9.23
## 12             9.37

Yanıt

Hamilelik sırasında sigara içen ve içmeyen annelerin bebeklerinde herhangi aşırı ya da beklenmedik doğum ağırlıkları var mıydı? Kısacası, bu sorunun cevabı evet. Sigara içme durumunu göz ardı edildiğinde, 31 doğum ağırlığı aykırı olarak belirlenmiştir. Doğum ağırlıklarının merkezi eğilimine kıyasla 4.75 poundun altındaki ve 10.25 poundun üstündeki değerler aşırıdır. Sigara içme durumu göz önünde bulundurulduğunda, hiç sigara içmemiş annelerin bebeklerinde doğum ağırlıkları daha tutarlıdır (daha az varyasyon). Sonuç olarak, yaklaşık 5.25 poundun altında ya da 10.00 poundun üzerinde olan bebekler aykırı olarak sınıflandırılmıştır. Hamilelik sırasında sigara içen annelerin bebeklerinin doğum ağırlıkları ortalama olarak daha düşüktür ama aynı zamanda daha tutarsızdır (daha fazla varyasyon). Bu nedenle, bu alt gruptaki az sayıda doğum ağırlığı aşırı olarak kabul edilir.

6. Mekansal Analiz

Soru

Amerika Birleşik Devletleri Nüfus Sayım Bürosu’nun (USCB) misyonu, “ülkenin insanları ve ekonomisi hakkında kaliteli veri sağlayan önde gelen kurum olarak hizmet vermek”tir. USCB web sitesi, Büro tarafından toplanan çeşitli verileri listeler, bunların bazıları Anayasa tarafından zorunlu kılınmıştır. İkamet eden nüfus bağlamında, aşağıdaki iki soru sorulduğunu varsayalım: Birleşik Devletler’in bitişik topraklarının nüfusu nasıl dağılmıştır? En yoğun nüfuslu eyaletin başkenti ile en az yoğun nüfuslu eyaletin başkenti arasındaki mesafe ne kadardır?

Veri

#load tidyverse and maps libraries
library(tidyverse)
library(maps)

## 
## Attaching package: 'maps'

## The following object is masked from 'package:purrr':
## 
##     map

#import USA map data
state_map <- map_data("state")

#plot USA map
ggplot(data=state_map,aes(x=long,y=lat,group=group)) +
  geom_polygon(color="black",fill="white") +
  labs(title="Lower 48 States of the USA",
       x="Longitude",
       y="Latitude",
       caption="Source: maps package") +
  theme_bw()

#load covidcast library
library(covidcast)

## We encourage COVIDcast API users to register on our mailing list:
## https://lists.andrew.cmu.edu/mailman/listinfo/delphi-covidcast-api
## We'll send announcements about new data sources, package updates,
## server maintenance, and new features.

#import census data
data(state_census)

#display data structure
glimpse(state_census)

## Rows: 57
## Columns: 9
## $ SUMLEV            <dbl> 10, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, …
## $ REGION            <chr> "0", "3", "4", "4", "3", "4", "4", "1", "3", "3", "3…
## $ DIVISION          <chr> "0", "6", "9", "8", "7", "9", "8", "1", "5", "5", "5…
## $ STATE             <dbl> 0, 1, 2, 4, 5, 6, 8, 9, 10, 11, 12, 13, 15, 16, 17, …
## $ NAME              <chr> "United States", "Alabama", "Alaska", "Arizona", "Ar…
## $ POPESTIMATE2019   <dbl> 328239523, 4903185, 731545, 7278717, 3017804, 395122…
## $ POPEST18PLUS2019  <int> 255200373, 3814879, 551562, 5638481, 2317649, 306175…
## $ PCNT_POPEST18PLUS <dbl> 77.7, 77.8, 75.4, 77.5, 76.8, 77.5, 78.1, 79.6, 79.1…
## $ ABBR              <chr> "US", "AL", "AK", "AZ", "AR", "CA", "CO", "CT", "DE"…

#import State data
data(state)

#create land area data frame
state_area <- data.frame(region=state.name,
                         area=state.area)

#display data structure
glimpse(state_area)

## Rows: 50
## Columns: 2
## $ region <chr> "Alabama", "Alaska", "Arizona", "Arkansas", "California", "Colo…
## $ area   <dbl> 51609, 589757, 113909, 53104, 158693, 104247, 5009, 2057, 58560…

#join and wrangle census and area data
state_density <- state_census %>%
  transmute(region=NAME,
            pop=POPESTIMATE2019) %>%
  inner_join(state_area,by="region") %>%
  transmute(region,density=pop/area) %>%
  filter(!region %in% c("Alaska","Hawaii"))

#display data structure
glimpse(state_density)

## Rows: 48
## Columns: 2
## $ region  <chr> "Alabama", "Arizona", "Arkansas", "California", "Colorado", "C…
## $ density <dbl> 95.00639, 63.89940, 56.82819, 248.98529, 55.24126, 711.77620, …

#change State name to lower case
state_density <- state_density %>%
  mutate(region=tolower(region))

#join density table to map table
state_popmap <- state_map %>%
  left_join(state_density,by="region")

Analiz

#create choropleth map of density
densmap <- ggplot() +
  geom_polygon(data=state_popmap,
               aes(x=long,y=lat,group=group,fill=density),
               color="black") +
  labs(title="Estimated Population Density of Contiguous US in 2019",
       fill="Residents\nper\nSquare\nMile",
       caption="Source: United States Census Bureau") +
  theme_void()

#display choropleth map
densmap

#load viridis library
library(viridis)

#display opposing color choropleth
densmap +
scale_fill_viridis(option="C")

#isolate highest and lowest density
state_density %>%
  arrange(density) %>%
  filter(row_number()==1 | row_number()== n())

##       region     density
## 1    wyoming    5.910891
## 2 new jersey 1133.510720

#create haversine function
haversine <- function(lat1,lon1,lat2,lon2){
  
  #convert to radians
  lat1 <- lat1*pi/180
  lat2 <- lat2*pi/180
  
  #convert differences to radians
  difflat <- (lat2-lat1)*pi/180
  difflon <- (lon2-lon1)*pi/180
  
  #compute subparts of formula
  sub1 <- (sin(difflat/2))^2
  sub2 <- (sin(difflon/2))^2
  sub3 <- sqrt(sub1+cos(lat1)*cos(lat2)*sub2)
  
  #compute and return distance
  dist <- 2*3958.8*asin(sub3)
  return(dist)
  
}

#compute haversine distance
haversine(41.16108,-104.80545,40.21705,-74.74294)

## [1] 1567.127

#compute haversine distance
haversine(47.60801,-122.33517,35.22709,-80.84312)

## [1] 2106.156

#create capital data frame
caps <- data.frame(city=c("Cheyenne","Trenton"),
                   lat=c(41.16108,40.21705),
                   long=c(-104.80545,-74.74294),
                   grp=c(1,1))

#display complementary choropleth
densmap +
scale_fill_viridis(option="C") +
  labs(subtitle="Distance between Capitals of Most and Least Dense States") +
  geom_point(data=caps,
             aes(x=long,y=lat),
             size=2,fill="red",shape=23) +
  geom_line(data=caps,
             aes(x=long,y=lat,group=grp),
            size=1,color="red",linetype="dashed") +
  geom_text(data=caps,
            aes(x=long,y=lat,label=city),
            nudge_x=-2,nudge_y=1,color="red") +
  annotate(geom="text",x=-92,y=39,label="1,567 miles",color="red")

Değerlendirme

#isolate lowest 10 densities
state_density %>%
  arrange(density) %>%
  head(n=10)

##          region   density
## 1       wyoming  5.910891
## 2       montana  7.263780
## 3  north dakota 10.784151
## 4  south dakota 11.482069
## 5    new mexico 17.234305
## 6         idaho 21.387376
## 7      nebraska 25.048338
## 8        nevada 27.864628
## 9        kansas 35.414203
## 10         utah 37.754463

Yanıt

ABD’nin bitişik eyaletlerindeki nüfus, kuzeydoğuda ve doğu ile batı kıyı şeridinde daha yoğun bir şekilde dağılmıştır. Kuzeydoğu eyaletleri arasında, New Jersey mil kare başına 1.134 sakin ile en yoğun nüfuslu eyalettir. Ortabatı ve merkez ABD’nin nüfus yoğunlukları, kıyı bölgelerine göre çok daha düşüktür. En düşük yoğunluklu eyaletlerin çoğu, Wyoming, Montana ve New Mexico gibi Rocky Dağları koridoru boyunca yer alır. Wyoming, mil kare başına sadece 6 sakin ile en düşüktür.

En yoğun nüfuslu başkent Trenton, New Jersey ile en az yoğun nüfuslu başkent Cheyenne, Wyoming arasındaki mesafe yaklaşık 1.567 mildir. Bu, belirli araştırma sorusuna cevap verirken, New Jersey ile karşılaştırıldığında benzer şekilde dramatik nüfus yoğunluğu azalmaları sunan Güney Dakota, Nebraska ve Kansas gibi daha yakın eyaletlerin de olduğu belirtilmelidir. Araştırmanın birincil amacı ABD’de daha az kalabalık bölgeleri keşfetmekse, bu daha yakın eyaletler yeterli olabilir.

7. Doğrusal İlişki

Soru

Palmer Station Uzun Süreli Ekolojik Araştırma (LTER) Programı, Antarktika’daki Palmer Takımadaları boyunca kutup deniz biyomunu inceler. Bu girişimlerden biri, yerel penguen popülasyonunun biyolojik özelliklerini incelemeyi içerir. Palmer pengueninin yüzgeç uzunluğu ile vücut kütlesi arasında bir ilişki var mı?

Veri

#load Palmer library
library(palmerpenguins)

#import penguin data
data(penguins)

#display structure of penguin data
glimpse(penguins)

## Rows: 344
## Columns: 8
## $ species           <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel…
## $ island            <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgerse…
## $ bill_length_mm    <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, …
## $ bill_depth_mm     <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, …
## $ flipper_length_mm <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186…
## $ body_mass_g       <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, …
## $ sex               <fct> male, female, female, NA, female, male, female, male…
## $ year              <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007…

#wrangle penguins data
penguins2 <- penguins %>%
  select(flipper_length_mm,body_mass_g) %>%
  na.omit()

#summarize penguin data
summary(penguins2)

##  flipper_length_mm  body_mass_g  
##  Min.   :172.0     Min.   :2700  
##  1st Qu.:190.0     1st Qu.:3550  
##  Median :197.0     Median :4050  
##  Mean   :200.9     Mean   :4202  
##  3rd Qu.:213.0     3rd Qu.:4750  
##  Max.   :231.0     Max.   :6300

Analiz

#create basic scatter plot
plot(penguins2$flipper_length_mm,
    penguins2$body_mass_g)

#create ggplot scatter plot
scatter <- ggplot(data=penguins2,
                  aes(x=flipper_length_mm,y=body_mass_g)) +
  geom_point() +
  labs(title="Palmer Archipelago (Antarctica)",
       subtitle="Penguin Measurements",
       x="Flipper Length (mm)",
       y="Body Mass (g)",
       caption="Source: Palmer Station LTER Program") +
  theme_bw()

#display scatter plot
scatter

#compute Pearson correlation coefficient
cor(penguins2$flipper_length_mm,
    penguins2$body_mass_g,
    method="pearson")

## [1] 0.8712018

#display best fit line
scatter +
  geom_smooth(method="lm",formula=y~x,se=FALSE) +
  geom_abline(slope=50,intercept=-4500,
              color="red",linetype="dashed") +
  geom_abline(slope=25,intercept=-1500,
              color="red",linetype="dashed")

#build linear model
linear_model <- lm(body_mass_g~flipper_length_mm,data=penguins2)

#extract linear model coefficients
coefficients(summary(linear_model))

##                      Estimate Std. Error   t value      Pr(>|t|)
## (Intercept)       -5780.83136 305.814504 -18.90306  5.587301e-55
## flipper_length_mm    49.68557   1.518404  32.72223 4.370681e-107

Değerlendirme

#display misleading best fit line
ggplot(data=penguins,
       aes(x=bill_depth_mm,y=body_mass_g)) +
  geom_point() +
  geom_smooth(method="lm",formula=y~x,se=FALSE) +
  labs(title="Palmer Archipelago (Antarctica)",
       subtitle="Penguin Measurements",
       x="Bill Depth (mm)",
       y="Body Mass (g)",
       caption="Source: Palmer Station LTER Program") +
  theme_bw()

#display multiple best fit lines
ggplot(data=penguins,
       aes(x=bill_depth_mm,y=body_mass_g,color=species,group=species)) +
  geom_point() +
  geom_smooth(method="lm",formula=y~x,se=FALSE) +
  labs(title="Palmer Archipelago (Antarctica)",
       subtitle="Penguin Measurements",
       x="Bill Depth (mm)",
       y="Body Mass (g)",
       color="Species",
       caption="Source: Palmer Station LTER Program") +
  theme_bw()

Yanıt

Palmer Takımadaları’nda 300’den fazla penguenin ölçümüne dayanarak, bir penguenin yüzgeç uzunluğu ile vücut kütlesi arasında güçlü, pozitif ve doğrusal bir ilişki olduğu görülmektedir. Özellikle, bir penguenin vücut kütlesinin, yüzgeç uzunluğundaki her bir milimetrelik artış için ortalama olarak yaklaşık 50 gram artması beklenmektedir. Bu ilişki tüm üç türde de tutarlıdır.

8. Lojistik İlişki

Soru

Murat Koklu ve arkadaşlarının yaptığı araştırma, dijital görüntülerle birlikte makine öğrenimi algoritmalarını kullanarak iki çeşit kabak çekirdeğini ayırt etmeyi amaçlamaktadır. Yazarlar, çekirdek tipini tahmin etmek için görüntülerden elde edilen çeşitli ölçümleri (örneğin, uzunluk, genişlik ve alan) uygulamaktadır. Sonuç olarak, modelleri, sadece ölçülen boyutlara dayanarak çekirdek tipini %87-89 doğrulukla ayırt edebilmektedir. Şöyle bir soru sorulduğunu varsayalım: Bir kabak çekirdeğinin uzunluğu ile türü arasında bir ilişki var mıdır?

Veri

#import pumpkin seed data
pumpkin <- read_csv("Pumpkin_Seeds_Dataset.csv")

## Rows: 2500 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): Type
## dbl (12): Area, Perimeter, Major_Axis_Length, Minor_Axis_Length, Convex_Area...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

#display structure of data
glimpse(pumpkin)

## Rows: 2,500
## Columns: 13
## $ Area              <dbl> 56276, 76631, 71623, 66458, 66107, 73191, 73338, 696…
## $ Perimeter         <dbl> 888.242, 1068.146, 1082.987, 992.051, 998.146, 1041.…
## $ Major_Axis_Length <dbl> 326.1485, 417.1932, 435.8328, 381.5638, 383.8883, 40…
## $ Minor_Axis_Length <dbl> 220.2388, 234.2289, 211.0457, 222.5322, 220.4545, 23…
## $ Convex_Area       <dbl> 56831, 77280, 72663, 67118, 67117, 73969, 73859, 704…
## $ Equiv_Diameter    <dbl> 267.6805, 312.3614, 301.9822, 290.8899, 290.1207, 30…
## $ Eccentricity      <dbl> 0.7376, 0.8275, 0.8749, 0.8123, 0.8187, 0.8215, 0.79…
## $ Solidity          <dbl> 0.9902, 0.9916, 0.9857, 0.9902, 0.9850, 0.9895, 0.99…
## $ Extent            <dbl> 0.7453, 0.7151, 0.7400, 0.7396, 0.6752, 0.7165, 0.71…
## $ Roundness         <dbl> 0.8963, 0.8440, 0.7674, 0.8486, 0.8338, 0.8480, 0.88…
## $ Aspect_Ration     <dbl> 1.4809, 1.7811, 2.0651, 1.7146, 1.7413, 1.7535, 1.64…
## $ Compactness       <dbl> 0.8207, 0.7487, 0.6929, 0.7624, 0.7557, 0.7522, 0.77…
## $ Type              <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A…

#wrangle pumpkin seed data
pumpkin2 <- pumpkin %>%
  transmute(Length=Major_Axis_Length,
            Type=as.factor(Type)) %>%
  na.omit()

#summarize pumpkin seed data
summary(pumpkin2)

##      Length      Type    
##  Min.   :320.8   A:1300  
##  1st Qu.:415.0   B:1200  
##  Median :449.5           
##  Mean   :456.6           
##  3rd Qu.:492.7           
##  Max.   :661.9

Analiz

#summarize stats by type
pumpkin2 %>%
  group_by(Type) %>%
  summarize(count=n(),
            average=mean(Length),
            stdev=sd(Length)) %>%
  ungroup()

## # A tibble: 2 × 4
##   Type  count average stdev
##   <fct> <int>   <dbl> <dbl>
## 1 A      1300    426.  37.2
## 2 B      1200    489.  54.9

#display box plot of length by type
ggplot(data=pumpkin2,aes(x=Length,y=reorder(Type,-Length))) +
  geom_boxplot(color="black",fill="sky blue") +
  labs(title="Pumpkin Seed Classification",
       subtitle="(based on digital imagery)",
       x="Major Axis Length (pixels)",
       y="Seed Type",
       caption="Source: Koklu et al. (2021)") +
  scale_x_continuous(limits=c(300,700),breaks=seq(300,700,50)) +
  theme_bw()

#create scatter plot of length versus type
bin_scatter <- ggplot(data=pumpkin2,aes(x=Length,y=ifelse(Type=="A",1,0))) +
  geom_jitter(width=0,height=0.05,alpha=0.2) +
  labs(title="Pumpkin Seed Classification",
       subtitle="(based on digital imagery)",
       x="Major Axis Length (pixels)",
       y="Proportion that are Type A",
       caption="Source: Koklu et al. (2021)") +
  scale_x_continuous(limits=c(300,700),breaks=seq(300,700,50)) +
  scale_y_continuous(limits=c(0,1),breaks=seq(0,1,0.1)) +
  theme_bw()

#display scatter plot
bin_scatter

#display line on scatter plot
bin_scatter +
  geom_smooth(method="lm",formula=y~x,se=FALSE)

#display sigmoid curve on scatter plot
bin_scatter +
  geom_smooth(method="glm",formula=y~x,se=FALSE,
              method.args=list(family="binomial"))

#build logistic model
logistic_mod <- glm(Type=="A"~Length,data=pumpkin2,
                    family="binomial")

#extract logistic model coefficients
coefficients(summary(logistic_mod))

##                Estimate  Std. Error   z value      Pr(>|z|)
## (Intercept) 13.43589960 0.561415083  23.93220 1.415952e-126
## Length      -0.02935021 0.001234749 -23.77019 6.794587e-125

#compute McFadden pseudo-R-squared
1-(summary(logistic_mod)$deviance/
   summary(logistic_mod)$null.deviance)

## [1] 0.2723096

Değerlendirme

Yanıt

2.500 kabak çekirdeğini analiz ettikten sonra, bir kabak çekirdeğinin ana eksen uzunluğu ile türü arasında orta derecede, negatif, doğrusal olmayan bir ilişki buluyoruz. Özellikle, Tip B çekirdekleri genellikle Tip A çekirdeklerinden daha uzundur. Sonuç olarak, bir çekirdeğin Tip A olma olasılığı, çekirdeğin uzunluğu arttıkça azalmaktadır.

9. Zaman Serisi Analizi

Soru

Boston Atletizm Birliği (BAA), 1887 yılında kurulmuş olup, 1897 yılından bu yana prestijli Boston Maratonu’nu düzenlemektedir. Her yıl binlerce koşucu, dünyanın en eski yıllık maratonunun şampiyonu olmak için yarışan dünya çapında profesyonel atletler de dahil olmak üzere, 26.2 mil uzunluğundaki parkuru tamamlar. Etkinliğin tarihine ilişkin daha fazla bilgi burada bulunabilir. Boston Maratonu’nu koşma süresi için hangi tarihsel eğilimler mevcuttur?

Veri

#import full html webpage
boston <- read_html("https://www.baa.org/races/boston-marathon/results/champions")

#extract tables from html webpage
boston_tabs <- html_nodes(boston,"table")

#isolate men's winning times
boston_men <- html_table(boston_tabs[[1]])

#display data structure
glimpse(boston_men)

## Rows: 126
## Columns: 4
## $ X1 <int> 2023, 2022, 2021, 2019, 2018, 2017, 2016, 2015, 2014, 2013, 2012, 2…
## $ X2 <chr> "Evans Chebet", "Evans Chebet", "Benson Kipruto", "Lawrence Cherono…
## $ X3 <chr> "Kenya", "Kenya", "Kenya", "Kenya", "Japan", "Kenya", "Ethiopia", "…
## $ X4 <chr> "2:05:54", "2:06:51", "2:09:51", "2:07:57", "2:15:58", "2:09:37", "…

#correct year 2023
boston_men$X1[1] <- 2023

#wrangle Boston data
boston_men2 <- boston_men %>%
  transmute(Year=as.integer(X1),
            Time=hms(X4))

#display data structure
glimpse(boston_men2)

## Rows: 126
## Columns: 2
## $ Year <int> 2023, 2022, 2021, 2019, 2018, 2017, 2016, 2015, 2014, 2013, 2012,…
## $ Time <Period> 2H 5M 54S, 2H 6M 51S, 2H 9M 51S, 2H 7M 57S, 2H 15M 58S, 2H 9M …

#compute minute per mile pace
boston_men3 <- boston_men2 %>%
  mutate(Miles=ifelse(Year<1924,24.5,26.2),
         Minutes=period_to_seconds(Time)/60,
         Pace=Minutes/Miles) %>%
  select(Year,Pace) %>%
  arrange(Year)

#display data structure
glimpse(boston_men3)

## Rows: 126
## Columns: 2
## $ Year <int> 1897, 1898, 1899, 1900, 1901, 1902, 1903, 1904, 1905, 1906, 1907,…
## $ Pace <dbl> 7.149660, 6.612245, 7.127891, 6.519728, 6.097279, 6.661224, 6.591…

Analiz

#display line graph
ggplot(data=boston_men3,aes(x=Year,y=Pace)) +
  geom_line(color="blue") +
  labs(title="Boston Marathon Winning Pace",
       subtitle="(no race in 2020)",
       x="Year",
       y="Winning Pace (min/mile)",
       caption="Source: Boston Athletic Association") +
  scale_x_continuous(limits=c(1897,2023),breaks=seq(1897,2023,14)) +
  scale_y_continuous(limits=c(4.5,7.25),breaks=seq(4.5,7.25,0.25)) +
  theme_bw()

#create ggplot scatter plot
ggplot(data=boston_men3,aes(x=lag(Pace,1),y=Pace)) +
  geom_point() +
  geom_smooth(method="lm",formula=y~x,se=FALSE) +
  labs(title="Boston Marathon Winning Pace",
       subtitle="(no race in 2020)",
       x="First Lag of Winning Pace (min/mile)",
       y="Winning Pace (min/mile)",
       caption="Source: Boston Athletic Association") +
  theme_bw()

#build lagged model
lag_model <- lm(Pace~lag(Pace,1),data=boston_men3)

#extract lagged model coefficients
coefficients(summary(lag_model))

##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)  0.7162600 0.21674443  3.304629 1.246207e-03
## lag(Pace, 1) 0.8661051 0.03928162 22.048608 1.490121e-44

#compute lagged model trend
boston_men4 <- boston_men3 %>%
  mutate(Trend=0.761626+0.8661051*lag(Pace,1))

#display trend on line graph
ggplot(data=boston_men4,aes(x=Year)) +
  geom_line(aes(y=Pace),color="blue") +
  geom_line(aes(y=Trend),color="red",linetype="dashed") +
  labs(title="Boston Marathon Winning Pace",
       subtitle="(blue=actual and red=predicted)",
       x="Year",
       y="Winning Pace (min/mile)",
       caption="Source: Boston Athletic Association") +
  scale_x_continuous(limits=c(1897,2023),breaks=seq(1897,2023,14)) +
  scale_y_continuous(limits=c(4.5,7.25),breaks=seq(4.5,7.25,0.25)) +
  theme_bw()

Değerlendirme

#display misleading line graph
ggplot(data=boston_men3,aes(x=Year,y=Pace)) +
  geom_line(color="blue") +
  labs(title="Boston Marathon Winning Pace",
       subtitle="(no race in 2020)",
       x="Year",
       y="Winning Pace (min/mile)",
       caption="Source: Boston Athletic Association") +
  scale_x_continuous(limits=c(1897,2023),breaks=seq(1897,2023,14)) +
  scale_y_continuous(limits=c(0,60),breaks=seq(0,60,10)) +
  theme_bw()

Yanıt

Boston Maratonu’nun 127 yıllık tarihinde, yarış kazananının ortalama koşu hızında belirgin bir azalma eğilimi görülmektedir. Sonuç olarak, yarış kazananlarının bitirme süreleri yıldan yıla düşmektedir. Bununla birlikte, dünyanın en iyi koşucularını çeken ödül parasının 1980’lerde mevcut hale gelmesinden sonra azalma hızı düzelmeye başlamıştır. O zamandan beri, bitirme süreleri insan performansının sınırlarına yaklaşmış gibi görünmekte ve büyük azalmalar çok daha az yaygın hale gelmiştir.

Keşfedici analize üç potansiyel genişletme hemen akla gelmektedir. İlk olarak, iki bilinen kırılma noktası (yarış mesafesi ve ödül parası) için daha dikkatli bir işlem yapılması muhtemelen gereklidir. Yarıştaki bu sistematik değişiklikleri hesaba katarak, kırılmalar arasındaki eğilimler için daha doğru modeller keşfedebiliriz. İkincisi, keşfi analizi çıkarımsal bir analize genişletebiliriz. Ayrıntılı nedenleri bu metnin kapsamı dışında olmakla birlikte, birinci dereceden otoregresyonun eğim parametresinin mutlak değerde 1’den az olması önemlidir. Eğim için en iyi tahminimiz \(\beta_1 = 0.866\)’dır, ancak tahmin için hata payını sağlamadık. Temel eğilim etrafındaki doğal değişkenlik göz önüne alındığında, gerçek eğimin 1’e (veya daha büyük) eşit olması mümkündür. 4. Bölümde, tahminler için hata payını hesaplamayı ve istatistiksel anlamlılıklarını değerlendirmeyi öğreniriz. Son bir genişletme olarak, ek gecikmeler veya diğer açıklayıcı değişkenler içeren daha karmaşık bir tahmini model geliştirebiliriz. Birden fazla değişken içerdikleri için bu tür modellere basit (tek değişkenli) doğrusal regresyon modellerine kıyasla çoklu doğrusal regresyon modelleri denir.